/*
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
#if !defined _gaussian_functor_
#define _gaussian_functor_

#include <src/algorithm/random/all.h>
//#include <src/math/Operators/GenericAlgebra.h>
//#include <src/miscellaneous/Function_Utilities.h>


/*template<class data, class result = data, std::size_t Dim = data::Dim>
struct Gaussian_Function {
	OS_STATIC_CHECK(false); // to implement
};*/

template<class Real_Type = OS_double>
struct Gaussian_Function 
	: public template_function_utility<OS_normal_distribution<Real_Type> > 
{
	typedef template_function_utility Base;
	Gaussian_Function(argument_type average = 0, argument_type sigma = 0.05) : Base(average, sigma)
	{};
};


/*! \brief multi-dimensional normal generator

	Base class (non compilable class) normal generator. 
*/
template<class data, class result = data, std::size_t Dim = data::Dim>
struct Gaussian_Functor {
	OS_STATIC_CHECK(false); // to implement
};

/*! \brief multi-dimensional normal generator

	We denote here \f$\sigma >0\f$ a real and consider the gaussian function of average \f$\overline{x}\f$ and variance \f$\sigma\f$
    \f[
        \mathcal{N}_{\sigma,\overline{x}}(x) = \frac{1}{\sqrt{2 \pi}\sigma} \exp\left(-\frac{(x-\overline{x})^2}{\sigma}\right)
    \f]
    for \f$x\in \mathcal{R}\f$.
	Defaults values are average \f$\overline{x} = 0\f$, \f$\sigma = 0.05\f$. 
*/
template<class data>
struct Gaussian_Functor<data,data, 2> 
{
	typedef typename data::value_type value_type;
	template <class I>
	Gaussian_Functor(const I& average, typename data::base_value_type sigma = 0.05) 
		: sigma_(sigma), average_(average.size())
	{
		std::copy(average.begin(),average.end(),average_.begin());
	};

	template <class bound>
	data operator()(const bound & nb_part){
		typedef boost::random::UniformEquallySpacedGenerator<value_type,2> uniform_generator;
		data temp(nb_part);
		uniform_generator local(nb_part);
		boost::random::GaussianGenerator<uniform_generator,2>	generator(average_, sigma_,local);
		generator.fill(temp);
		return temp;
	};
	value_type average_;
	typename data::base_value_type sigma_;
};


template<class data>
struct Gaussian_Functor<data,data,1> 
{
	Gaussian_Functor(typename data::value_type average = 0.5, typename data::base_value_type sigma = 0.05) 
		: average_(average),sigma_(sigma) {};
	data operator()(OS_size nb_part){
		typedef boost::random::UniformEquallySpacedGenerator<data::value_type> uniform_generator;
		data temp(nb_part);
		uniform_generator local(nb_part);
		boost::random::GaussianGenerator<uniform_generator,1>	generator(average_, sigma_,local);
		generator.fill(temp);
		return temp;
	};
	typename data::value_type average_;
	typename data::base_value_type sigma_;
};

template <class geometric_process, OS_size Dim = geometric_process::Dim >
struct geometric_process_transport_map {
	typedef geometric_process									stochastic_process_type;
	typedef typename stochastic_process_type::underlying_type	underlying_type;
	typedef typename stochastic_process_type::matrix_type		matrix_type;
	OS_STATIC_CHECK(false); // not yet implemented
};

/*! Optimal Transport Map for the geometric brownian motion.
/* Optimal Transport Map for the geometric brownian motion.
	This is the one dimensional case.
*/
template <class geometric_process>
struct geometric_process_transport_map<geometric_process,0> 
{
	typedef geometric_process stochastic_process_type;
	typedef typename stochastic_process_type::underlying_type	underlying_type;
	typedef typename stochastic_process_type::matrix_type		matrix_type;
	geometric_process_transport_map() 
	{
	};
	geometric_process_transport_map(const smart_ptr<stochastic_process_type> &process, OS_double time) : time_(time)
	{
		x0_ = process->stochastic_process_type::OS_Base::initial_conditions();
		drift_ = process->stochastic_process_type::OS_Base::drift(time,x0_); // this is (r_+1/2 \sigma ^2)
		matrix_type		volat = process->stochastic_process_type::OS_Base::volatility(time,x0_);
		InverseCumulativeGaussian_ = OS_inverse_cumulative_normal_distribution<underlying_type>(0.,volat*std::sqrt(time));
	};
	underlying_type operator()(const underlying_type & uniform_sampling) 
	{
		return x0_*std::exp(InverseCumulativeGaussian_(uniform_sampling)+drift_*time_);
	};
	underlying_type x0_;
	underlying_type drift_; // this is (r_1/2 \sigma ^2)
	OS_inverse_cumulative_normal_distribution<underlying_type> InverseCumulativeGaussian_;
	OS_double time_;
};

/*! Optimal Transport Map for the Cox Ingelsson Ross process.
/* Optimal Transport Map for the Cox Ingelsson Ross process.
	This class throw a compile error : CIR processes are not defined for higher dimensions.
*/
template <class CIR_process, OS_size Dim = CIR_process::Dim >
struct CIR_process_transport_map {
	typedef CIR_process											stochastic_process_type;
	typedef typename stochastic_process_type::underlying_type	underlying_type;
	typedef typename stochastic_process_type::matrix_type		matrix_type;
	OS_STATIC_CHECK(false); // not yet implemented
};

/*! Optimal Transport Map for the Cox Ingelsson Ross process.
/* Optimal Transport Map for the Cox Ingelsson Ross process.
	This is the one dimensional case.
*/
template <class CIR_process>
struct CIR_process_transport_map<CIR_process,0> : public operator_base<
																		typename CIR_process::underlying_type,
																		typename CIR_process::underlying_type
																	>
{
	typedef CIR_process											stochastic_process_type;
	typedef typename stochastic_process_type::underlying_type	underlying_type;
	typedef typename stochastic_process_type::matrix_type		matrix_type;
	CIR_process_transport_map() 
	{
	};
	CIR_process_transport_map(const smart_ptr<stochastic_process_type> &process, OS_double time) : time_(time)
	{
			OS_double kappa = process->kappa_;
			OS_double x0	= process->S0_;
			OS_double sigma	= process->sigma_;
			OS_double theta = process->theta_;
			delta_			= process->delta_;
			mu_				= (4.*kappa)/(sigma*sigma*(1.-exp(-kappa*time_)));
			lambda_			= mu_*exp(-kappa*time_)*x0;
	};
	underlying_type operator()(const underlying_type & uniform_sampling) const
	{
			OS_inverse_cumulative_non_central_chi_squared_distribution<OS_double> chi(delta_,lambda_);
			return chi(uniform_sampling)/mu_;
	};
	OS_double time_;
	OS_double mu_;
	OS_double delta_;
	OS_double lambda_;
};


template<class data, std::size_t = data::Dim>
struct InverseCumulativeFunctor : public operator_base<data,data>
{
	data operator()(OS_size nb_part, OS_size dimension){
		data temp(nb_part,dimension);
		typedef CompoundGenerator<boost::random::sobol_generator<32>, OS_inverse_cumulative_normal_distribution<> > generator_type;
		generator_list<data::value_type,generator_type > generator(dimension);
		generator.fill(temp);
//		std::cout << std::endl << "---InverseCumulativeFunctor---" << temp << std::endl;
		return temp;
	};
};

/*! \brief one dimensional inverse normal cumulative generator

    Let \f$\sigma >0\f$ a real and consider the gaussian function of average \f$\overline{x}\f$ and variance \f$\sigma\f$
    \f[
        \mathcal{N}_{\sigma,\overline{x}}(x) = \frac{1}{\sqrt{2 \pi \sigma}} \exp\left(-\frac{(x-\overline{x})^2}{\sigma}\right)
    \f]
    for \f$x\in \mathcal{R}\f$. Let \f$s^{-1}(x)=\int_{-\infty}^x \mathcal{N}_{\sigma}(u)du\f$ be the normal cumulative and denote \f$s(y)\f$ its inverse.

    This generator return a data vector \f$\{v_i \in \mathcal{R}\}_{i=1,..,N}\f$ filled with
    \f[
        v_i = s(y_i)
    \f]
    where \f$\{ y_i \in [0,1]\}_{i=1,..,N}\f$ is described by an uniform generator :
    \f[
        y_i = \frac{i+.5}{N+1}
    \f]
    Defaults values are average \f$\overline{x} =0\f$, variance \f$\sigma = 0.05\f$.
*/
template<class data> 
struct InverseCumulativeFunctor<data,1>
{
	InverseCumulativeFunctor(OS_double a = 0. , OS_double s = std::sqrt( 0.1 )) : function_(a,s){};
	data operator()(OS_size nb_part){

		typedef boost::random::UniformEquallySpacedGenerator<OS_double,0> uniform_generator_type;
		uniform_generator_type generator_(nb_part);
		data temp(nb_part);
		generator_.fill(temp);
		return function_.operator ()<data,data>(temp);
	};
	OS_inverse_cumulative_normal_distribution<typename data::base_value_type> function_;
};






#endif